function[temp] = rob_regress1(X,Y,Z,X1,Y1),
     [rs,rstat] = robustfit(X,Y);
     [rs1,rstat1] = robustfit(X,Z);
     [rs2,rstat2] = robustfit(X,X1);
     residuals = abs(rstat.resid);
     residuals1 = abs(rstat1.resid);
     residuals2 = abs(rstat2.resid);
     [y,i] = sort(residuals);
     [y1,i1] = sort(residuals1);
     [y2,i2] = sort(residuals2);
     
     rej_fac = floor(length(X)*0.2);%The last value is the percent of data to be rejected
%      temp = i(end-20:end);
     temp = intersect(i(end-rej_fac:end),i1(end-rej_fac:end));
     temp = intersect(temp,i2(end-rej_fac:end));
     % The number 20 is ar
    % Robust outlier rejection  --- END